function [UCoef,ExtraInfo]...
        =Bellman_U(PP,SS,Pir,T,Profit,TAU,ir,w,UCoefp,Flag_SS)
% Note:
%   This Bellman equation maps from EVFunCoefp to EVFunCoef.
% Profit=SS.Profit;ir_dom=SS.ir_dom;ir_ext=SS.ir_ext;Pir=SS.Pir;dE=SS.dE;
% T=SS.T; w_N=SS.w_N; w_H=SS.w_H; State=SS.FunApp.V.State;TAU=SS.TAU;
%% Sequential Solver
% Step 1: from UCoefp to VbarCoef
[Opt_Vbar,Vbar]     =   Fun_Vbar(PP,SS,T,Profit,TAU,ir,w,UCoefp,SS.FunApp.V.State,Flag_SS);
VbarCoef_Mat        =   SS.FunApp.V.BasMat_0\reshape(Opt_Vbar,SS.FunApp.V.MatSize_State);
VbarCoef            =   VbarCoef_Mat(:);
% Step 2: from VCoef to VbarCoef
[Opt_U,U]           =   Fun_U(PP,SS,Pir,VbarCoef,SS.FunApp.U.State,Flag_SS);
UCoef_Mat           =   SS.FunApp.U.BasMat_0\reshape(Opt_U,SS.FunApp.U.MatSize_State);
UCoef               =   UCoef_Mat(:);
%% Collocation Operation
AccType     =   'Collocation';
if Flag_SS==1 
    switch AccType
        case 'Collocation'
            % Collocation
            % VBasMat*VbarCoef = FlowVec+CoefMat*UCoefp;
            FlowVec     =   Vbar.FlowVec;
            CoefMat     =   Vbar.CoefMat;
            % from V to U
            FlowVec     =   U.FlowVec+U.CoefMat*( SS.FunApp.V.BasMat_0_State\FlowVec );
            CoefMat     =   U.CoefMat*( SS.FunApp.V.BasMat_0_State\CoefMat );
            
            UCoef       =   ( SS.FunApp.U.BasMat_0_State-CoefMat )\FlowVec;
        case 'Howard'
            % Howard
            HowardItNum     =   40;
            for ii=1:HowardItNum
                Temp_V      =   Vbar.FlowVec+Vbar.CoefMat*UCoefp;
                Temp_VbarCoef=  SS.FunApp.V.BasMat_0_State\Temp_V;
                Temp_U      =   U.FlowVec+U.CoefMat*Temp_VbarCoef;
                UCoef       =   SS.FunApp.U.BasMat_0_State\Temp_U;

                Diff        =   max(abs(UCoef-UCoefp));
                if Diff<1e-8
                    break;
                end
                UCoefp     =   UCoef;
            end
    end
end
%% Collect the Information
ExtraInfo   =   struct('FunCoef',struct('Vbar',VbarCoef,'U',UCoef),...
                       'Vbar',Vbar,'U',U);
